/*
  Copyright (C) 2011, Kyungjoo Kim
  
  This file is part of LINAL (LINear ALgebra)
  
  All rights reserved.

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  
  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.
  
  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  Also add information on how to contact you by electronic and paper mail.
  
  
  Kyungjoo Kim
  iamkyungjoo@gmail.com
*/

#include "environment.hxx"

int main(int argc, char **argv) {

  linal::initialize( argc, argv );

  bool help;
  int nthreads, datatype, m, n, n_rhs, b_mn, locker, balancer;
  double tol;

  linal::get_option_bool  ( "-help", "list options",
                            __FILE__,LINAL_OPTION_MODE_RUN, help);

  linal::get_option_int   ( "-datatype", LINAL_OPTION_TEXT_DATATYPE,
                            __FILE__,TEST_DATATYPE, datatype, help);

  linal::get_option_int   ( "-M", "size of matrix",
                            __FILE__,N, m, help);
  linal::get_option_int   ( "-N", "size of matrix",
                            __FILE__,N, n, help);
  linal::get_option_int   ( "-Nrhs", "# of rhs",
                            __FILE__,1, n_rhs, help);

  linal::get_option_int   ( "-b_mn", "block size of matrix",
                            __FILE__,B_MN, b_mn, help);

  linal::get_option_int   ( "-locker", "locker size",
                            __FILE__,24, locker, help);

  linal::get_option_int   ( "-balancer", "balancer",
                            __FILE__,0, balancer, help);

  linal::get_option_int   ( "-nthreads", "# of threads",
                            __FILE__,1, nthreads, help);

  linal::get_option_double( "-tol", "error tolerence",
                            __FILE__,LINAL_ERROR_TOL, tol, help);

  if (help == LINAL_OPTION_MODE_DRY) {
    linal::finalize();
    return 0;
  }

  omp_set_num_threads(nthreads);

  linal::host::Manager host = new linal::host::Manager_();
  host->initialize( balancer );

  linal::Dispatcher dispatcher = new linal::Dispatcher_();
  linal::Dispatcher_::Group_ group(b_mn, locker);

  dispatcher->add_device( host );
  dispatcher->add_group( group );

  dispatcher->lock();
  std::cout << *dispatcher << std::endl;

  linal::set_dispatcher( dispatcher );

  linal::Flat_ A[2], B, X[2], p[2], norm;
  linal::Hier_ HA, hp;

  for (int i=0;i<2;++i) {
    A[i].create( datatype, m, n );
    p[i].create( FLA_INT, m, 1 );
    X[i].create( datatype, m, n_rhs );
  }

  B.create( datatype, n, n_rhs );
  norm.create( linal::get_datatype_real(datatype), 1, 1 );

  HA.create( A[0], b_mn, b_mn );
  hp.create( p[0], b_mn, b_mn );

  FLA_Random_matrix( ~A[0] );
  FLA_Copy( ~A[0], ~A[1] );

  FLA_Random_matrix( ~B );

  if (linal::get_verbose()) {
    FLA_Obj_show( "= A[0] =", ~A[0], "%4.3f", "===" );
    FLA_Obj_show( "= A[1] =", ~A[1], "%4.3f", "===" );
  }

#pragma omp parallel
  {
#pragma omp single nowait
    linal::lu_incpiv( HA, hp );
  }

  FLA_LU_piv( ~A[1],~p[1] );

  for (int i=0;i<2;++i)
    FLA_LU_piv_solve( ~A[i], ~p[i], ~B, ~X[i] );

  if (linal::get_verbose()) {
    FLA_Obj_show( "= A[0] =", ~A[0], "%4.3f", "===" );
    FLA_Obj_show( "= p[0] =", ~p[0], "%3d", "===" );

    FLA_Obj_show( "= A[1] =", ~A[1], "%4.3f", "===" );
    FLA_Obj_show( "= p[1] =", ~p[1], "%3d", "===" );

    FLA_Obj_show( "= X[0] =", ~X[0], "%4.3f", "===" );
    FLA_Obj_show( "= X[1] =", ~X[1], "%4.3f", "===" );
  }

  FLA_Axpy( FLA_MINUS_ONE, ~X[0], ~X[1] );
  FLA_Norm1( ~X[1], ~norm );

  if (linal::get_verbose())
    FLA_Obj_show( "=Diff X[1]=", ~X[1], "%4.3f", "===" );

  float  norm_f = *((float *)FLA_Obj_buffer_at_view( ~norm ));
  double norm_d = *((double*)FLA_Obj_buffer_at_view( ~norm ));

  printf("- TEST::");
  for (int i=0;i<argc;++i)
    printf(" %s ", argv[i]);
  printf("\n");

  int rval;
  switch (linal::get_datatype_real(datatype)) {
  case LINAL_SINGLE_REAL:
    printf("SINGLE PRECISION\n");
    if (norm_f < LINAL_ERROR_TOL) {
      printf("PASS::Norm :: %10.9f \n", norm_f);   rval = 0;
    } else {
      printf("FAIL::Norm :: %10.9f \n", norm_f);   rval = -1;
    }
    break;
  case LINAL_DOUBLE_REAL:
    printf("DOUBLE PRECISION\n");
    if (norm_d < LINAL_ERROR_TOL) {
      printf("PASS::Norm :: %E \n", norm_d);   rval = 0;
    } else {
      printf("FAIL::Norm :: %E \n", norm_d);   rval = -1;
    }
    break;
  }

  // ---------------------------------------
  // ** Matrices  
  HA.free();
  hp.free();

  for (int i=0;i<2;++i) {
    A[i].free();
    p[i].free();
    X[i].free();
  }

  B.free();
  norm.free();

  linal::finalize();

  return rval;
}
